program novel
use mod_part_interp
implicit none
integer,parameter :: n1=32,n2=32,Np=1,Nstokes=2
double precision,parameter::cfl=0.01d0,tmax=1000.
double precision :: pi
double precision :: dx,dy,x,y,allocated_memory,dt,time
double precision :: length,uamp
double precision :: x1,y1,vx1,vy1,uxanalytic=0.,uyanalytic=0.
double precision,dimension(n1,n2) :: ux,uy
integer :: i,j
logical :: lgradu=.false.
character*100 :: cnpar,formp
!--------------------------
pi=4.*datan(1.0d0)
length=2.*pi
uamp=0.1
dx=length/dfloat(n1)
dy=length/dfloat(n2)
open(unit=1,file='vel.out',status='unknown')
do j=1,n2
  y = dy*dfloat(j-1)
  do i=1,n1
    x = dx*dfloat(i-1)
    ux(i,j) = uamp*sin(y) 
    uy(i,j) = uamp*sin(x)
!    write(*,*)i,j,x,y
!    write(1,*)x,y,ux(i,j),uy(i,j)
  enddo
enddo
close(1)
!------------
call allocate_particles(Np,Nstokes,lgradu,allocated_memory)
call init_particle(n1,n2,dx,dy,ux,uy)
dt=(cfl)*(dx/uamp)
time=0.
write(*,*) 'starting time marching ...'
formp = '('//trim(adjustl('4'))//'es16.3E1)'
do while(time.lt.tmax)
!  write(*,*) 't=',time
!  call calc_particle_diagnostic()
  call get_particle_coordinate(1,1,x1,y1,vx1,vy1)
  write(*,*)time,x1,y1,vx1,vy1
!  uxanalytic=uamp*sin(y1)
!  uyanalytic=uamp*sin(x1)
!  write(*,*)'analytic:',uxanalytic,uyanalytic
  call evolve_particles_euler(ux,uy,n1,n2,dt)
  time=time+dt
enddo
write(*,*) '.........done'
!------------
endprogram novel
